Method to generate the in-situ state of stress in a domain ω in a geological structure

ABSTRACT

The invention relates to a method and a computer-implemented invention for numerical modeling of a geological structure. The present invention solves the problem providing a method for use in the numerical simulation of the in-situ stress in a geological structure represented by a domain Ω located under its external ground surface S. The method comprises mainly two steps: determining a first state of in-situ stress in the domain Ω by means of six stress components and a second step determining a correction of the first state of stress in order to satisfy the equilibrium equation.

OBJECT OF THE INVENTION

The invention relates to a method and a computer-implemented invention for numerical modeling of a geological structure.

PRIOR ART

Geological materials located at a certain depth are subjected to stresses, whether from gravitational or tectonic origin or residual stresses from previous geological activity.

Gravitational stress is induced by the weight of the overburden; the residual stress is contributed by strain energy stored in the rock during the geological past; and, the tectonic stress is caused by current tectonic forces affecting the region. When measurements are made at a particular site, the instruments sense the combined contribution of these three components.

Greenfield conditions are the un-modified initial stresses that exist within the ground, arising from the various geological processes that it has been subjected to. These stresses are modified by construction or any other man-made activity.

When an opening is excavated in these materials, or the conditions are changed in some other way, e.g. pore pressure variations induced by fluid injection or extraction in a reservoir, the local stress field is modified so that a new set of stresses appear in the rock mass.

Knowledge of the magnitudes and directions of both pre-existing (greenfield) and man-induced stresses is an essential component of any geomechanical analysis, since for instance if the existing stresses values are near the limit strength of the material, the new stress may more easily reach that strength and potentially lead to failure or collapse.

Together with the details of geologic structure, in-situ stresses existing in the ground constitute one of the most important starting data for the analysis. Developing stresses depend essentially upon in-situ stresses, and on pore pressures also related to the in-situ situation. A good definition of these in-situ stresses is therefore essential.

Following the generation of the finite element mesh, and before any future activities can be analyzed, the initial conditions in the ground must therefore be established.

In the most general context, this could be ideally achieved by modeling the complete geological history of the site, if that is known. However, if that is not known, or the necessary analysis is too complicated, the initial conditions may be achieved by various different simplified procedures, which may lead to significantly different in-situ stress states. Because of its important consequences for the subsequent analysis, it is therefore important to understand and select the right procedure for initialization of in-situ stresses in Finite Elements calculations.

Throughout the document, the stress tensor is expressed as a stress vector, because the computer methods that usually carry out the invention manipulate this stress tensor in this vector form. However, the method according to the invention can be carried out with other ways of expressing the stress tensor, as the method relates to the operations that are realized over its components.

The invention presented herein is a new such method of estimating the in-situ stress conditions with no need of simulating the evolution of the geological history.

DESCRIPTION OF THE INVENTION

The present invention solves the problem providing a method for use in the generation of an estimate of the in-situ stress field in a geological structure represented by a domain Ω. The method comprises the steps:

-   a) determining a first state of in-situ stress, representable as the     vector σ^(prop) expressible by means of the stress components     σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz), in the domain Ω;     -   This first step a) starts with an estimated guess of the in-situ         stress field that satisfies the desired vertical/horizontal         stress ratio or any other condition that is considered to affect         the in-situ stresses. The stress ratios may be obtained from         tectonic stress data bases and even corrected with additional         local data. When determining a first state of in-situ stress by         means of the stress components it must be interpreted that the         method is formulated in the more general scenario, i.e. a         three-dimensional domain. If for instance a two dimensional         domain is used then only existing variables will be involved in         the calculation; that is, (σ_(xx), σ_(zz), σ_(xz)). The other         horizontal stress σ_(yy) exists but it can be calculated from         existing variables as a result of a postprocessing step. The         same interpretation applies to a coordinate such as x=(x₀, y₀,         z₀) being in the two-dimensional case modeled only as x=(x₀,         z₀).     -   The estimated guess obtained as a result of step a) is not         always in equilibrium. Therefore, the method comprises a         subsequent step b) providing an equilibration stage by imposing         forces to the numerical mesh that are equal and opposite to the         nodal forces that reflect the stress imbalance. The unbalanced         nodal forces (residual forces) are then applied and         redistributed in a single calculation (in the case of a linear         elastic material), or according to an embodiment, in an         iterative process until its components are sufficiently small         (in the case of a non-linear elastic material). This step b) is         carried out by means of sub-stages b.1) to b.9): -   b) determining a correction σ^(corr) of the first state of in-situ     stress σ^(prop) carrying out the following steps:     -   b.1) determining a finite element discretization of the domain Ω         by means of a numerical mesh, being Ω_(e) the domain of a         particular finite element;     -   b.2) determining the boundary conditions on the boundary B of         the numerical mesh of the domain Ω, comprised of Dirichlet or         prescribed displacements on a part B₁ of the boundary B, and a         Newmann or prescribed stresses on the remaining part B₂ of the         boundary, being B=B₁ ∪ B₂;     -   b.3) determining the gravity equivalent element force f_(e)         ^(ext) in the numerical mesh for each element using the         expression

f_(e) ^(ext)=∫_(Ω) _(e) N^(T) ρg dV

-   -   -   wherein N is the element shape function matrix, ρ is the             mass density of the material and g is the gravity vector;

    -   b.4) determining the internal element forces f_(e) ^(int)         corresponding to the first state of in-situ stress obtained in         the preceding step a) in the finite element mesh for all         elements from the expression

f_(e) ^(int)=∫_(Ω) _(e) B^(T) σ^(prop) dV

-   -   -   wherein B is the element matrix relating strains and nodal             displacements, and σ^(prop) is the stress tensor expressible             as (σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz))^(T);

    -   b.5) determining the element stress imbalance f_(e) ^(res) as         the subtraction between the gravity equivalent element force and         the internal element forces

f _(e) ^(res) =f _(e) ^(ext) −f _(e) ^(int)

-   -   b.6) assembling the element residual forces f_(e) ^(res)         representing the stress imbalance at the element level into a         global residual forces vector f^(res), representing the stress         imbalance in each node of the numerical mesh at the domain         structural level; that is, each component of the global residual         forces vector f^(res) comprises the sum, extended over all         elements, of the components of element residual forces f_(e)         ^(res) related to the same node of said component of the global         residual forces vector f^(res); and, loading the domain         discretized by the numerical mesh with the residual forces         vector f^(res) and solving the global system of equations

Ku=f^(res)

-   -   -   where K is the global stiffness matrix and u is the global             nodal displacements vector, with K resulting from the             assembly of the individual elements stiffness matrix K_(e)             that corresponds to the integral over the volume of each             finite element in the mesh:

K_(e)=∫_(Ω) _(e) B^(T) EB dV

-   -   -   and the assembly procedure consists of the sum on the global             stiffness matrix of the contribution of each individual             stiffness matrix considering the transformation between the             local numbering of each element and the global numeration at             the structural level that accounts for the correspondence of             degrees of freedom as it is standard procedure in finite             element formulation. In the previous expression E is the             elasticity stiffness matrix with components E_(ijkl) that             define the linear elastic behavior σ_(ij)=E_(ijkl)Δε_(kl) ,             wherein σ_(ij) is the stress tensor, ε_(kl) is the strain             tensor and E_(ijkl) is the elasticity tensor,

    -   b.7) determining the deformation by means of the compatibility         equation as

$ɛ_{kl} = {\frac{1}{2}\left( {u_{k,l} + u_{l,k}} \right)}$

-   -   b.8) determining the correction of the first state of in-situ         stress by means of the material equation

Δσ_(ij) ^(corr)=E_(ijkl)Δε_(kl)

-   -   -   where E_(ijkl) is the same tensor of elasticity at each             point of the domain, which has been used for the calculation             of the global stiffness matrix K in step b.6).

    -   b.9) provide the sought in-situ stress σ as the correction of         the first state of in-situ stress given in step a)         σ^(prop)=(σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz))         as=σ^(prop)+Δσ^(corr).

    -   After the correction, an in-situ stress field in equilibrium is         obtained.

DESCRIPTION OF THE DRAWINGS

These and other features and advantages of the invention will be seen more clearly from the following detailed description of a preferred embodiment provided only by way of illustrative and non-limiting example in reference to the attached drawings.

FIG. 1 This figure shows a scheme with the main steps of the method according to the invention.

FIG. 2 This figure shows the domain Ω used in a first embodiment to apply the method of the invention.

FIG. 3 This figure shows the Boundary Conditions used in a first embodiment of the method of the invention.

FIG. 4 This figure shows the in-situ stress field which results of the application of the step a) of a method of the invention in a first embodiment.

FIG. 5 This figure shows the in-situ stress field which results of the application of the step b) of a method of the invention in a first embodiment.

FIG. 6 This figure shows the principal in-situ stress components along a vertical section of the domain in FIG. 2.

FIG. 7 This figure shows the domain Ω used in a first and a second embodiments to apply the method of the invention (particular examples 2 and 3)

FIG. 8 This figure shows the in-situ stress field which results of the application of the step a) of a method of the invention in a first embodiment (particular example 2)

FIG. 9 This figure shows the in-situ stress field which results of the application of the step b) of a method of the invention in a first embodiment (particular example 2)

FIGS. 10a-10d These figures show the in-situ stress components along four vertical cross-sections of the domain in FIG. 7 in a first embodiment (particular example 2).

FIG. 11 This figure shows the in-situ stress field which results of the application of the step a) of a method of the invention in a second embodiment (particular example 3).

FIG. 12 This figure shows the in-situ stress field which results of the application of the step b) of a method of the invention in a second embodiment (particular example 3).

FIG. 13a-13d These figures show the in-situ stress components along four vertical cross-sections of the domain in FIG. 7 in a second embodiment (particular example 3).

DETAILED DESCRIPTION OF THE INVENTION

General Approach

In general terms, the set of governing equations commonly used for the mathematical solution of the stress-strain problem in a rock mass is composed of: (a) the Equilibrium Equations, (b) the Compatibility Equations, and (c) the Constitutive or Material Laws. In Cartesian coordinates and under the assumption of small strains, these laws are:

-   -   a) Equilibrium: σ_(ij,j)+ρ_(i)=0     -   b) Compatibility: ε_(kl)=½(u_(k,l)+u_(l,k))     -   c) Material law: σ_(ij)=f(ε_(kl), ξ) in general, or         σ_(ij)=E_(ijkl)ε_(kl) for linear elastic materials, wherein is         one or more internal variables (also called history variables)         that may be present in the specific formulation used for         inelastic material laws.

The previous set of partial differential equations, which are valid over the domain Ω limited by the boundary B, must be complemented with appropriate boundary conditions, in a general case in the form of:

-   -   d) Dirichlet or prescribed u_(i)=u_(i) on the part B₁ of the         boundary where the displacements are known, and     -   e) Newmann σ_(ij)n_(j)=t_(i) on the remaining part B₂ of the         boundary (B=B₁ ∪ B₂) where the external forces are known.

In the case of in-situ stress calculations, which includes strictly speaking only the evaluation of the stress state, in the most relevant cases the only equations that must be strictly satisfied in the whole domain Ω are the equilibrium equations (a), together with the corresponding Newmann condition (e) on the part of the boundary B₂. The Material law (c) on its side must be satisfied sometimes and partially; only when the material is a non-linear elastic material with a limit stress condition (i.e. strength criterion or plastic stress limit), and in that case neither the strength criterion nor the plastic stress limit can be violated by the in-situ stress solution obtained.

In spite of this, the evaluation of the in-situ stress state in the prior art normally involves the solution of the whole set of equations (a)-(e) with some real or fictitious material parameters and some Dirichlet Boundary Conditions, although the resulting displacements and strains are normally discarded (and reset to zero values) after the calculation. The reasons for solving the whole system are:

-   -   The system of equilibrium equations (a)+Newmann Boundary         Conditions (e), with or without the constraints given by the         stress limits of the material law (c), does not constitute a         complete set of equations. That is, there are more unknowns than         equations, and therefore such system has an infinite number of         solutions that correspond to the many possible in-situ stress         states that may exist on a given domain, depending for instance         on geological and tectonic history.     -   Finite Elements codes are normally only prepared to solve the         full system of equations (a)-(e), that is with certain values of         material parameters and Dirichlet Boundary Conditions on B₁. The         solutions obtained in this way for stresses indeed satisfy         equilibrium conditions (a) and (e).     -   The particular material laws and parameters used for the in-situ         stress calculations do not have to correspond necessarily to a         realistic material behavior. They may be assigned fictitious         values, which are chosen for convenience in order to generate         the desired effects, such as for instance the desired ratio of         vertical to horizontal stresses K₀=σ_(h)/σ_(v).

Concerning the validity of the assumption of small strain, it is also true that in its previous geological history until the current in-situ state, the material may have been subject to large strains. However, the usual methods for the calculation of in-situ stress do not try to follow the stress-strain history of the rock material through its geologic history (in general very complicated and often uncertain), but to obtain the current stress state on the basis of an estimate which satisfies mathematically the equilibrium conditions and, on the other hand, approaches as much as possible the remaining information that is available on the in-situ stress distribution.

As for pore pressures, in a particular embodiment in-situ pore pressures are evaluated separately based on relatively simple assumptions, generally hydrostatic distributions, with reference values based on available information. Some particular examples are those wherein the pore pressure is obtained from advanced models based on field measures or seismic data.

In some embodiments where that is required for the in-situ stress Finite Elements calculation (i.e. if the material is assumed as a non-linear elastic material) those pore pressures are then introduced as fixed prescribed values at each point, in the iterative calculation process between strain and (total) stress.

Simplified procedures similar to the ones used for the pore pressures may be employed for parts of the domain that may represent geological materials that exhibit viscous incompressible flow in the long term, such as salt domes, since that behavior leads to stress fields that are also practically hydrostatic.

A particular embodiment of a method according to the invention comprises the following steps:

Step a) Initial Estimate of the In-Situ Stress State

The process starts providing an estimate σ^(prop) of the existing in-situ stress state, expressible by means of the stress components σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz). This does not need to be an exact evaluation, and the stress field proposed in this step does not need to strictly verify the equilibrium conditions of equation (a), since a correction step will be performed later. However, it is convenient that the estimate:

-   -   satisfies as exactly as possible the characteristics of the         in-situ field as known from field measurements or other sources         (e.g. World-Stress-Map, www.world-stress-map.com);     -   is as close as possible to equilibrium, so that the variations         that will take place after the correction step will be         minimized, and therefore the final equilibrated state will         deviate as little as possible from the proposed stress state

To unambiguously assess the closeness of this estimate to equilibrium, the lower the residual given by the expression μσ_(ij,j) ^(prop)+ρg_(i)∥ is, the closer is the estimate to the equilibrium.

According to this invention, three alternative embodiments of the method are considered to generate the first estimate of in-situ stresses:

-   -   1) Taking the vertical stress as a function of the depth of the         corresponding element or Gauss point, with simultaneous estimate         of the horizontal stress ratio based on an input value of         parameter K₀. Therefore, in this case the estimate of the         in-situ stress will exhibit vertical and horizontal principal         directions, and the following stress values:         -   a. σ_(xx)=σ_(H), σ_(yy)=σ_(h), σ_(zz)=σ_(V), σ_(xy)=0,             σ_(xz)=0, σ_(yz)=0         -   b. where σ_(V)=∫γ dz, σ_(h)=K₀σ_(v) and σ_(H)=K_(an)σ_(h),             wherein γ is the specific weight of the material, K₀ is the             ratio between the minimum horizontal stress and vertical             stress, and K_(an) is a predetermined horizontal stress             ratio estimating the horizontal stress anisotropy.     -   In case of two-dimensional analysis, the stress values are, in a         simplified form, the following:         -   a. σ_(zz)=∫γdz, σ_(xx)=K₀σ_(zz), τ_(xz)=0         -   b. where γ is the specific weight of the material, and K₀             the desired horizontal-to-vertical stress ratio.     -   Former integral values σ_(v)=∫γdz, or σ_(zz)=∫γdz may be         generalized when the vertical stress is known at certain depth         σ_(v) ₀ wherein the calculation are expressed as         σ_(v)=∫γdz+σ_(v) ₀ , or σ_(zz)=∫γdz+σ_(v) ₀ being the integral         extended along a path defined from the depth where σ_(v) ₀ is         taken until the depth where the σ_(V) or σ_(zz) is calculated.     -   2) Subdividing geometrically the domain of interest in a         collection of subdomains, and using for each of them a         minimization procedure involving the so-called stress functions,         which are scalar potential functions leading to stress fields as         their derivatives.         -   For general three-dimensional analysis, these scalar             potential functions are the so-called Maxwell stress             functions A(x,y,z), B(x,y,z), C(x,y,z) from which the             various stress components follow:

$\mspace{20mu} {\sigma_{xx} = {\frac{\partial^{2}B}{\partial z^{2}} + \frac{\partial^{2}C}{\partial y^{2}} + {K_{0\; x}\gamma \; z} + \sigma_{{xx}\; 0}}}$ $\mspace{20mu} {\sigma_{yy} = {\frac{\partial^{2}C}{\partial x^{2}} + \frac{\partial^{2}A}{\partial z^{2}} + {K_{0\; y}\gamma \; z} + \sigma_{{yy}\; 0}}}$ $\mspace{20mu} {\sigma_{zz} = {\frac{\partial^{2}A}{\partial y^{2}} + \frac{\partial^{2}B}{\partial x^{2}} + {\gamma \; z} + \sigma_{{zz}\; 0}}}$ ${\sigma_{xy} = {{- \frac{\partial^{2}C}{{\partial x}{\partial y}}} + \sigma_{{xy}\; 0}}};{\sigma_{xz} = {{- \frac{\partial^{2}B}{{\partial x}{\partial z}}} + \sigma_{{xz}\; 0}}};{\sigma_{yz} = {{- \frac{\partial^{2}A}{{\partial y}{\partial z}}} + \sigma_{{yz}\; 0}}}$

-   -   -   wherein K_(0x) and K_(0y) are the predefined ratios of             horizontal-to-vertical stresses in the x- and y-directions             respectively, σ_(xx0), σ_(yy0), σ_(zz0), σ_(xy0), σ_(xz0)             and σ_(yz0) represent the values of the in-situ stress state             at the point x=y=z=0 of the domain Ω, and wherein A(x,y,z),             B(x,y,z), C(x,y,z) are three functions that must obey the             additional following condition:

${{\nabla^{4}A} + {\nabla^{4}B} + {\nabla^{4}C}} = {\frac{3}{2 - v}\left( {\frac{\partial^{2}A}{\partial x^{2}} + \frac{\partial^{2}B}{\partial y^{2}} + \frac{\partial^{2}C}{\partial z^{2}}} \right)}$

-   -   -   under the boundary conditions, wherein v is the Poisson             ratio.         -   In the particular case of a two-dimensional analysis, the             stress components to be specified and explicitly involved in             the solution of the partial differential equations are only             three (σ_(xx), σ_(zz), σ_(xz)). In this case the scalar             potential function is known as the Airy stress function             φ(x,z) from which the stresses are derived:

${\sigma_{xx} = {\frac{\partial^{2}\phi}{\partial z^{2}} + {K_{0}\gamma \; z} + \sigma_{{xx}\; 0}}};{\sigma_{zz} = {\frac{\partial^{2}\phi}{\partial x^{2}} + {\gamma \; z} + \sigma_{{zz}\; 0}}};{\sigma_{xz} = {{- \frac{\partial^{2}\phi}{{\partial x}{\partial z}}} + \sigma_{{xz}\; 0}}}$

-   -   -   wherein φ(x,z) is a function determined from the boundary             conditions, which automatically satisfies the differential             equations of equilibrium within the domain. (Note that if             φ(x, z) is selected equal to a constant or a first order             polynomial, the above expressions revert into the previous             estimate of method 1).

    -   In both cases (2D and 3D), a third degree polynomial form is         considered for the stress functions, with constant coefficients         A_(i) in 2D, and A_(i), B_(j), and C_(k) in 3D. By applying the         previous second derivative expressions these coefficients lead         to linear functions of the geometrical coordinates x,z (in 2D         case) or x,y,z (in 3D case) for the various components of the         stress tensor. These linear expressions involve some of the         coefficients A_(i), B_(j) or C_(k) of the polynomial stress         functions, which need to be determined and are the unknowns of         the subsequent minimization process.

    -   In full 3D, according to an embodiment, the stress functions are         selected as the following third degree polynomial:

${A\left( {x,y,z} \right)} = {{\frac{A_{1}}{6}x^{3}} + {\frac{A_{2}}{2}x^{2}y} + {\frac{A_{3}}{2}x^{2}{z++}\frac{A_{4}}{6}y^{3}} + {\frac{A_{5}}{2}y^{2}x} + {\frac{A_{6}}{2}y^{2}z} + {\frac{A_{7}}{6}z^{3}} + {\frac{A_{8}}{2}z^{2}x} + {\frac{A_{9}}{2}z^{2}y}}$ ${B\left( {x,y,z} \right)} = {{\frac{B_{1}}{6}x^{3}} + {\frac{B_{2}}{2}x^{2}y} + {\frac{B_{3}}{2}x^{2}{z++}\frac{B_{4}}{6}y^{3}} + {\frac{B_{5}}{2}y^{2}x} + {\frac{B_{6}}{2}y^{2}z} + {\frac{B_{7}}{6}z^{3}} + {\frac{B_{8}}{2}z^{2}x} + {\frac{B_{9}}{2}z^{2}y}}$ ${C\left( {x,y,z} \right)} = {{\frac{C_{1}}{6}x^{3}} + {\frac{C_{2}}{2}x^{2}y} + {\frac{C_{3}}{2}x^{2}{z++}\frac{C_{4}}{6}y^{3}} + {\frac{C_{5}}{2}y^{2}x} + {\frac{C_{6}}{2}y^{2}z} + {\frac{C_{7}}{6}z^{3}} + {\frac{C_{8}}{2}z^{2}x} + {\frac{C_{9}}{2}z^{2}y}}$

-   -   which, substituted in the previous expressions of the stress         components, lead to:

σ_(xx)=(B ₈ +C ₅)x+(B ₉ +C ₄)y+(B ₇ +C ₆ −K _(0x)γ)z+σ _(xx0)

σ_(yy)=(C ₁ +A ₈)x+(C ₂ +A ₉)y+(C ₃ +A ₇ −K _(0y)γ)z+σ _(yy0)

σ_(zz)=(A ₅ +B ₁)x+(A ₄ +B ₂)y+(A ₆ +B ₃−γ)z+σ _(zz0)

σ_(xy) =−C ₂ x−C ₅ y+σ _(xy0)

σ_(xz) =−B ₃ x−B ₈ z+σ _(xz0)

σ_(xy) =−A ₆ y−A ₉ z+σ _(xy0)

-   -   In 2D, according to a further embodiment the general format         considered for the Airy function is a third degree polynomial,

${\phi \left( {x,z} \right)} = {{\frac{A_{1}}{6}x^{3}} + {\frac{A_{2}}{6}x^{2}z} + {\frac{A_{3}}{6}{xz}^{2}} + {\frac{A_{4}}{6}z^{3}}}$

-   -   and the resulting stresses components of first degree in x and         z:

σ_(zz) =A ₁ x+(A ₂γ)z+σ _(zz0)

σ_(xx) =A ₃ x+(A ₄ +K ₀γ)z+σ _(xx0)

σ_(xz) =−A ₂ x−A ₃ z+σ _(xz0)

-   -   Some of the unknown coefficients (A₁, A₂, A₃, A₄ in 2D, or         A_(i), B_(j) and C_(k), with i, j, k=1,9 in 3D) can be         determined on the basis of fundamental considerations such as         the fact that the vertical gradients of vertical stress be γ and         of horizontal stress be K_(0y)γ (A₂=A₄=0 in 2D, and in 3D         B₇=C₆=A₇=C₃=0, A₆=−B₃) and B₉ +C₄, C₁+A₈, A₅ +B₁ and A₄ +B₂ can         be considered a single variable since they appear together in         the equations. The remaining unknown parameters are grouped in a         vector denoted as X_(i). (1=1,2 in 2D, and i=1,9 in 3D).     -   In a further embodiment, in order to determine these remaining         unknowns, a prismatic subdomain is considered that is limited on         the top and bottom by surfaces S₁ and S₂ which are plane but not         necessarily horizontal, and are subject to the following         boundary conditions:         -   Normal stress on the top surface is prescribed as a linear             function of x, y and z         -   Shear stress intensity on the top surface is also prescribed             as a linear function of x, y and z.         -   Shear stress on the bottom surface is linked to the amount             of normal stress on the same surface.     -   On the basis of those conditions, the following objective         function φ(X_(i)) is defined:

${\Phi \left( X_{i} \right)} = {{\int_{S_{1}}{\left( {\sigma_{n}^{(\alpha)} - {\overset{\_}{\sigma}}^{(\alpha)}} \right)^{2}\ {S_{1}}}} + {\int_{S_{1}}{\left( {\tau^{(\alpha)} - {\overset{\_}{\tau}}^{(\alpha)}} \right)^{2}\ {S_{1}}}} + {\int_{S_{2}}{\left( {\tau^{(\beta)} - {\overset{\_}{\tau}}^{(\beta)}} \right)^{2}\ {S_{2}}}}}$

-   -   The function is the minimized with respect to each of the         coefficients:

$\frac{\partial{\Phi \left( X_{i} \right)}}{\partial X_{i}} = {{{\frac{\partial}{\partial X_{i}}{\int_{S_{1}}{\left( {\sigma_{n}^{(\alpha)} - {\overset{\_}{\sigma}}^{(\alpha)}} \right)^{2}\ {{S_{1}++}}\frac{\partial}{\partial X_{i}}{\int_{S_{1}}{\left( {\tau^{(\alpha)} - {\overset{\_}{\tau}}^{(\alpha)}} \right)^{2}\ {S_{1}}}}}}} + {\frac{\partial}{\partial X_{i}}{\int_{S_{2}}{\left( {\tau^{(\beta)} - {\overset{\_}{\tau}}^{(\beta)}} \right)^{2}\ {S_{2}}}}}} = 0}$

-   -   In these expressions, σ_(n) ^((α)), τ^((α)) are the normal and         shear stresses of the proposed distribution on the top surface         of the subdomain, S₁; τ^((β)) is shear stress on the bottom         surface of the subdomain, S₂; σ ^((α)), τ ^((α)) are the normal         and shear stress values to be prescribed on S₁, and τ ^((β)) is         the shear stress value to be prescribed on S₂.

The above system of partial derivatives leads to the same number of equations as unknowns X_(i), which constitutes the algebraic system to be solved for the subdomain In order to obtain the X_(i).

The above process is sequentially applied to each of the subdomains from the surface to the bottom, applying for the top layer of each, the same conditions for normal stress as obtained at the common surface of the upper subdomain, and the same shear values as prescribed to the upper subdomain on the common surface.

-   -   3) A third procedure is also proposed, involving first the         elastic solution of the entire domain under gravity loads, and         then a recalculation of the horizontal stress component         according to the desired K₀. The detailed steps are the         following:         -   (i) to solve the elastic problem with gravity loads and on             the entire domain of interest         -   (ii) from this solution, to keep the values of the vertical             stress at each point of the domain, and         -   (iii) to recalculate horizontal stress values at each point             by means of the desired K₀ and K_(ani) coefficients.

Step b) Equilibrating the Proposed State

With the procedures disclosed in former step a), the estimated in-situ stress field will satisfy the desired stress ratio, but in general will not be in equilibrium. To ensure the consistency of the analysis, the method according to the invention includes an automatic equilibration stage b) which comprises imposing forces to the numerical mesh and the discretization given by steps b.1) and b.2) that are equal and opposite to the nodal forces that reflect the stress imbalance.

The unbalanced nodal forces (residual forces) are applied and redistributed in a single calculation (in the case of linear elastic materials) or with an iterative process until its components are sufficiently small (in the case of non-linear elastic materials, e.g. elasto-plastic).

The calculation of the unbalanced forces is performed as follows:

-   -   b.3) determining the gravity equivalent element force f_(e)         ^(ext) in the finite element numerical mesh for each element         using the expression

f_(e) ^(ext)=∫_(Ω) _(e) N^(T) ρg dV

-   -   -   wherein N is the element shape function matrix, ρ is the             mass density of the material and g is the gravity vector;

    -   b.4) determining the internal element forces f_(e) ^(int)         corresponding to the previously estimated state of stress from         the expression

f_(e) ^(int)=∫_(Ω) _(e) B^(T) σ^(prop) dV

-   -   -   wherein B is the element matrix relating strains and nodal             displacements, and σ^(prop) is the stress tensor expressible             as (σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz))^(T);

    -   b.5) determining the element stress imbalance (or residual         forces) f_(e) ^(res) as the subtraction between the gravity         equivalent element force and the internal element forces

f _(e) ^(res) =f _(e) ^(ext) −f _(e) ^(int)

-   -   b.6) assembling the element residual forces into the global         nodal residual forces vector f^(res), imposing the residual         forces as nodal loads acting on the finite element mesh and         solving the global system of equations

Ku=f^(res)

-   -   -   wherein K is the global stiffness matrix and u is the global             nodal displacements vector; and,

    -   next steps comprise accumulating the results of the previous         operation to the previously estimated state of in-situ stress;         that is,

    -   b.7) determining the deformation by means of the compatibility         equation as

ε_(kl)=½(u _(k,l) +u _(i,k))

-   -   b.8) determining the correction of the first state of in-situ         stress by means of the material equation

Δσ_(ij) ^(corr)=E_(ijkl)Δε_(kl)

-   -   -   where E_(ijkl) is the same tensor of elasticity at each             point of the domain, which has been used for the calculation             of the global stiffness matrix K in step b.6); and,

    -   b.9) providing the sought in-situ stress σ as the correction of         the first state of in-situ stress given in step a)         σ^(prop)=(σ_(xx), σ_(yy), σ_(zz), σ_(xy), σ_(xz), σ_(yz)) as         σ=σ^(prop)+Δσ^(corr).

In this manner, an in-situ stress field in equilibrium that approaches the initial guess of in-situ stress (e.g. imposed K₀ condition) is obtained.

In view of FIG. 1, a first estimate of in-situ stress field (100) is obtained by:

-   -   generating a guess for the in-situ stress field (101);     -   performing computations to ensure that the in-situ stress field         is in equilibrium (102); and,     -   accepting the computed equilibrated in-situ stress field as         valid and proceed (103) with the subsequent steps of a numerical         simulation based on the previous in-situ stress field result of         performing computations to ensure that the initial stress field         is in equilibrium (102).

The calculation of the unbalanced nodal forces and the correction of the initial estimate of in-situ stress are carried out by an iterative process if the material is a non-linear elastic material, e.g. elasto-plastic material, until the components of unbalanced forces, the corrections calculated from the unbalanced forces or both are sufficiently small.

PARTICULAR EXAMPLE 1

FIGS. 2 to 6 illustrate a first particular embodiment, wherein the domain Ω comprises a plurality of layers depth D and length 6D. FIG. 2 represents just the half of the domain Ω, since it is symmetric. The layers are in-depth and an overburden pressure σ_(ov)=γD is taken into account.

In the more general case, the vertical stress σ_(V) according to the direction z of gravity is determined according to the following steps:

-   -   determining the straight vertical path C connecting the         coordinate x and the point of the vertical projection of said         coordinate x at the upper surface S limiting the domain Ω,     -   determining σ_(V)=∫_(C) γ(z) dz wherein γ(z) is the specific         weight of the material and the integral sign denotes the line         integral along the path C,

If the vertical stress σ_(V) ₀ is known at certain depth, the path C starts on that depth and the integral is expressed as ∫_(C) γ(z) dz+σ_(V) ₀ . Expressed in an alternative manner, when an overburden exists on the upper surface S of the domain Ω, the boundary conditions comprises said overburden pressure; therefore, for determining σ_(V)=∫_(C) γ(z) dz along the path C the integral is evaluated at least in two parts, a first part for the contribution of the overburden at the coordinate x at the surface S in the path C; and a second part for the rest of the path C.

According to the invention, the upper surface S is the upper surface of the domain Ω. In most cases, the upper surface is the external ground surface. If this is not the case, the soil located above the upper surface S can be deemed as a mass generating an overburden pressure over the upper surface S.

The boundary conditions for this example are shown in FIG. 3, considering the vertical stress given by overweight, including overburden, and the horizontal stress given by desired values, different for each layer.

The method of the invention is carried out according to the described steps. The result of applying step a) is the in-situ stress field shown in FIG. 4. Step b) is carried out over this in-situ stress field to achieve the estimate shown in FIG. 5. The final in-situ stress field does not perturb the vertical stress, and does not dilute the desired horizontal stresses that are specified on the boundary.

In all these figures, the results are represented as two lines in each Gauss point representing the principal stress direction (line orientation) and modulus (line length). The vertical one represents the value of the vertical stress and the horizontal one represents the value of the horizontal stress.

FIG. 6 shows the principal in-situ stress components (vertical and horizontal stress) along a vertical section of the domain depicted in FIG. 2.

PARTICULAR EXAMPLE 2

FIG. 7 represents a second particular embodiment, wherein the domain Ω comprises a plurality of layers covering a depth D and length 3D. The domain of this particular example is not symmetric, since the interfaces between the layers are inclined, not parallel to the horizontal boundaries of the domain.

The proposed in-situ stress state is estimated by applying the step a) of the method according to the procedure 1) described above, that is, the values of parameters γ and K₀ for each layer which are shown in FIG. 7, are used directly to obtain the proposed stress a).

The result of step a) is the proposed in-situ stress field shown in FIG. 8. After this estimate, step b) of the method is applied, obtaining the in-situ stress field shown in FIG. 9.

FIGS. 10a, 10b, 10c and 10d show the components xx, zz and xz of the resulting in-situ stress along four vertical cross-sections of the domain depicted in FIG. 7. The sections represented are at 5, 45, 135 and 175 meters from the left boundary of the domain. The results show that the desired horizontal to vertical stress ratio K₀ is effectively recovered at all vertical cross-sections of the domain.

This in-situ stress, which turns out generally oriented with the geological layers, shows however some changes in direction near the vertical limits at left and right ends of the domain.

PARTICULAR EXAMPLE 3

In this example, the domain Ω is the same as that in the previous example, but the trial estimate of the in-situ stresses is obtained by means of the Airy function of the embodiment 2) of the step a) as shown in FIG. 11. Due to the additional constants, not only values of parameters γ and K₀ are enforced, but also the condition of no shear stress on the inter-layer contact surface.

The final in-situ stress field after applying the step b) of a method according to the invention to this estimate of the in-situ stress field can be observed in FIG. 12 and the stress components σ_(xx), σ_(zz), σ_(xz) are plotted in vertical cross sections in FIGS. 13a, 13b, 13c and 13d . The sections represented are at 5, 45, 135 and 175 meters from the left boundary of the domain.

As it can be seen in FIG. 12, the stresses also turn out generally oriented with the layers, but in this case this happens in the entire domain including the left and right boundaries of the domain. 

1. A method implemented in a computer for the numerical simulation of a porous medium that is capable of comprising hydraulic fracture, wherein the whole behavior of the porous medium is simulated following the following steps: defining a domain; creating a numerical mesh on the domain to generate a geomechanical and fluid model of the porous medium with standard finite elements comprising delimitation associated with the domain wherein edges of the mesh follow at least the boundaries of the domain and the pre-existing fractures, populating the numerical mesh elements with mechanical properties, at least elastic properties; and flow properties, at least permeability properties, the said populated numerical mesh forming a geomechanical and fluid flow model, populating the numerical mesh elements with the initial stress and initial fluid pressure values, assigning boundary conditions, assembling a system of equations, solving an evolution in time for a coupled solution of the fluid-flow model and the geomechanical models; wherein creating the numerical mesh comprises: creating a numerical mesh using zero-thickness interface elements which are inserted along a path of the pre-existing fractures, locating surfaces extending along finite element boundaries along which new fractures may propagate; and, creating a numerical mesh using zero-thickness interface elements along the mid-plane of the new fractures.
 2. The method according to claim 1 wherein the method comprises a constitutive geomechanical model for the zero-thickness elements located in the fractures, the constitutive geomechanical model being formulated in terms of: a total stress vector σ_(mp)=(σ_(n), τ_(l1), τ_(l2)) at the facture mid-plane where (n, l1, l2) is an orthogonal reference system with l1 and l2 aligned with the fracture mid-plane and n its normal; a displacement vector û=(u_(n), v_(l1), v_(l2)) defined as the relative displacement of two faced points in the zero-thickness elements located in the fracture, a hyperbolic failure surface F(σ_(n), τ) being τ=√{square root over (τ_(l1) ²+τ_(l2) ²)} defining the yield surface; wherein the F=0 condition results in a hyperbolic curve in the σ−τ plane, F=τ²−tan² φ(σ²−2ασ)=0 having two asymptotes intersecting with angle φ in respect to the σ axis in point α for τ=0; wherein given the τ_(l1) ^(n) and τ_(l1) ^(n) values in a predetermined time step n, the τ_(l1) ^(n+1) and τ_(l2) ^(n+1) values in the next time step n+1 are calculated as satisfying the following equation: ${\tan \left( \frac{\beta - \lambda}{2} \right)} = {{\tan \left( \frac{\beta - \lambda_{0}}{2} \right)}\left( \frac{\tau + \sqrt{\tau^{2} + {a^{2}\tan^{2}\varphi}}}{\tau_{0} + \sqrt{\tau_{0}^{2} + {a^{2}\tan^{2}\varphi}}} \right)^{- \frac{K_{t}\mspace{31mu} \Delta \; v}{K_{n}\Delta \; u_{n\mspace{14mu}}\tan \; \varphi}}}$ given that τ²−tan ²φ(σ²−2ασ)=0, wherein β relates imposed tangential relative displacements Δv_(l1) and Δv_(l2), K_(n) and K_(t) are normal and shear stiffness moduli respectively, λ₀ is the ratio between τ_(l1) ^(n) and τ_(l1) ^(n); and, λ is the ratio between τ_(l1) ^(n+1) and τ_(l1) ^(n+1).
 3. The method according to claim 1, wherein after solving the evolution in time for a coupled solution of the fluid-flow model and the geomechanical, the numerical mesh is remeshed according to displacements of the fractures to model the propagation of said fractures.
 4. The method according to claim 1, wherein the domain comprises at least an injection borehole and wherein before assembling the system of equations the method comprises assigning loads according to injection conditions.
 5. The method according to claim 1, wherein governing equations for the porous medium are: Equilibrium equation: L _(c) ^(T) σ+ρg=0 Effective stress principle: σ′=σ+α_(c) m _(c) p Constitutive relationship: σ′=D_(c)ε Compatibility equation: ε=L_(c)u where σ=(σ₁₁, σ₂₂, σ₃₃, σ₁₂, σ₁₃, σ₂₃) is a 6-component total stress vector, $L_{c} = \begin{pmatrix} \frac{\partial}{\partial x_{1}} & 0 & 0 & 0 & \frac{\partial}{\partial x_{3}} & \frac{\partial}{\partial x_{2}} \\ 0 & \frac{\partial}{\partial x_{2}} & 0 & \frac{\partial}{\partial x_{3}} & 0 & \frac{\partial}{\partial x_{1}} \\ 0 & 0 & \frac{\partial}{\partial x_{3}} & \frac{\partial}{\partial x_{2}} & \frac{\partial}{\partial x_{1}} & 0 \end{pmatrix}^{T}$ is a differential operator, g is a gravity vector, ρ is a bulk or averaged density of the fluid-solid system, σ′=(σ′₁₁, σ′₂₂, σ′₃₃, σ′₁₂, σ′₁₃, σ′₂₃) is a 6-component effective stress vector, α_(c) is a Biot-Willis coefficient, m_(c) is an equivalent to an identity tensor in this vector formulation, p is a pore fluid pressure, D_(c) is a rock constitutive matrix according to elastic conditions, ε is a strain vector, and u is a displacements vector; governing equations for the fluid flow in the continuum medium are: Continuity equation: ${{{div}\mspace{14mu} q} + {\left( {\frac{\alpha_{c} - \varnothing}{K_{s}} + \frac{\varnothing}{K_{f}}} \right)\frac{\partial p}{\partial t}} + {\alpha_{c}\frac{\partial ɛ_{v}}{\partial t}}} = 0$ Darcy's law (general form): $q = {{- \frac{K_{c}}{\gamma_{f}}}\left( {{{grad}\mspace{14mu} p} - {\rho_{f}g}} \right)}$ where q is a fluid flux, K_(s) is a solid bulk modulus, K_(f) is a fluid bulk modulus, Ø is porosity, ε_(v) is a volumetric strain, K_(c) is a rock hydraulic conductivity tensor, γ_(f) is a fluid specific weight; and ρ_(f) is a fluid density, equations governing mechanical behavior of fractures are posed at the fracture mid-plane and in a reference system aligned with the fracture mid-plane: Effective stress principle: σ′_(mp)=σ_(mp)+α_(f) m _(mp) p _(mp) Constitutive relationship: dσ′_(mp)=D_(mp)dû Compatibility equation: dû=L_(mp)du where σ′_(mp) is an effective stress vector at the fracture mid-plane, α_(f) is a fracture Biot's coefficient, m_(mp) is fluid pressure in a direction normal to a discontinuity axis, p_(mp) is a fluid pressure at the mid-plane of the fracture, û is a relative displacement vector, and L_(mp) is an operator that relates the relative displacement vector to an actual displacement at the faces of the fracture.
 6. The method according to claim 5, wherein the equations governing the fluid flow along the fracture are posed at the fracture mid-plane and in a reference system aligned with the fracture mid-plane: Continuity equation: ${{{div}\mspace{14mu} Q_{1}} + {\frac{1}{M_{f}}\frac{\partial p_{mp}}{\partial t}} + {\alpha_{f}\frac{\partial u_{n}}{\partial t}}} = 0$ Darcy's law: $Q_{1} = \left. {{- T_{1}}\mspace{11mu} {grad}} \middle| {}_{f}\left( {\frac{p_{mp}}{\gamma_{f}} + z} \right) \right.$ Transversal flow: q _(t) =K _(t)(p ^(bot) −p ^(top)) where u_(n) is the fracture aperture, Q₁ is a local flow rate vector along the fracture mid-plane, α_(f) and M_(f) respectively are a fracture Biot's coefficient and modulus, T_(l) is a mid-plane longitudinal transmissivity matrix, q_(t) is a local flow rate vector transversal to the fracture mid-plane, K_(t) is a transversal conductivity and p^(bot) and p^(top) are fluid pressures on each face of the fracture.
 7. The method according to claim 5, wherein the rock hydraulic conductivity is defined according to the volumetric strain (e_(v)) as: $K_{c} = {\left( {\frac{K_{i}}{n_{i}^{3}}\frac{\left( {^{\varepsilon_{v}} - \left( {1 - n_{i}} \right)} \right)^{3}}{^{\varepsilon_{v}}}} \right)I}$ where K_(i) is an initial hydraulic conductivity, n_(i) is an initial porosity, and I an identity matrix.
 8. The method according to claim 6, wherein the transmissivity along the fracture mid-plane is a cubic law where T_(l) is proportional to the cube of the aperture u_(n) ³.
 9. The method according to claim 1, wherein a region of the domain is populated with surfaces discretized using zero-thickness interface elements extending along the finite elements boundaries.
 10. The method according to claim 9, wherein the region of the domain populated with surfaces discretized using zero-thickness interface elements extending along the finite elements boundaries is the whole domain.
 11. A data processing system comprising means configured to carrying out a method according to claim
 1. 12. A computer program product adapted to perform a method according to claim
 1. 13. The method according to claim 4, wherein injection conditions include a flow-rate or fluid pressure. 